
estimate.mu <- function(data, sigma.gold, mean.mu = 0, sd.mu = 1, plot = T) {

	n = length(data)
	avg = sum(data)/n

	mu = seq(avg-3, avg+3, by=0.01)
	
	log.likelihood = rep(0, length(mu))
	mu.mle = rep(0, n)
	
	log.priori = -(mu - mean.mu)^2 * sigma.gold^2 / sd.mu^2
	log.posteriori = log.priori
	mu.map = rep(0,n)

	for(i in 1:n) {
		x = data[i]

		# MLE
		log.likelihood = log.likelihood - (mu - x)^2
		mu.mle.cur = mu[which.max(log.likelihood)]
		mu.mle[i] = mu.mle.cur
	
		# MAP
		log.posteriori = log.posteriori - (mu - x)^2	
		mu.map.cur = mu[which.max(log.posteriori)]
		mu.map[i] = mu.map.cur

		# plot
		if (plot == T) {
			par(mfrow=c(1,2))
			plot(mu,log.likelihood,type='l', ylab = "obj. func", xlab = "mu",
				main=paste('x[', i, '] = ', x,'\n','mu.mle = ', mu.mle.cur, '\n',
											'mu.map = ', mu.map.cur,sep=''))
			lines(mu,log.posteriori,col='blue')

			plot(1:i, mu.mle[1:i],type='l',
				xlim=c(1,n),ylim=c(avg-sigma.gold/2,avg+sigma.gold/2),
				xlab='sample size', ylab='mu.est')
			lines(1:i, mu.map[1:i],col='blue')

			readline()
			#Sys.sleep(0.001)
		}
	}
	return(list(mu.mle = mu.mle, mu.map = mu.map))
}

#n = 1000
#mu.gold = 3
#sigma.gold = 10
#mean.mu = 2.5
#sd.mu = 1

#mse.mle = rep(0,n); mse.map = rep(0,n)

#m = 100
#for (i in 1:m) {
#	data = rnorm(n, mean=mu.gold, sd = sigma.gold)
#	mu.est = estimate.mu(data, sigma.gold, mean.mu = mean.mu, sd.mu = sd.mu, plot=F)
#	mse.mle = mse.mle + (mu.est$mu.mle - mu.gold)^2
#	mse.map = mse.map + (mu.est$mu.map - mu.gold)^2
#}

#mse.mle = sqrt(mse.mle) / m
#mse.map = sqrt(mse.map) / m
#plot(1:n, mse.mle, type='l',ylim=c(min(min(mse.mle),min(mse.map)),max(max(mse.mle),max(mse.map))),
#	xlab = 'sample size', ylab = 'MSE')
#lines(1:n, mse.map, col='blue')


